©1995,1996, Richard B. Winston
Please Register
Last modified on Wednesday, January 31, 1996 - 7:12:26 PM
I recommend Anderson and Woessner (1992) to anyone who wishes to use groundwater models. It deals with many aspects of groundwater modeling not covered in this document.
The basic steps in creating a groundwater model are as follows:
- Develop a conceptual model of the system.
- Design a grid to represent the model area with as much detail as you need.
- Assign hydrologic properties to each grid cell.
- Select an appropriate time step for transient models.
- Set up boundary conditions.
- Select a solver.
- Run the model and compare with measured data to evaluate the effectiveness of the model.
- If necessary revise the model.
If groundwater modeling is part of some larger project it is a good idea to start the modeling effort at an early stage. Even a poorly calibrated model can provide insight into the system being modeled and may help in identifying deficiencies in the data which is collected. It would be disheartening, for example, to collect data for a year only to realize at the end that crucial data was missing that could easily have been collected earlier in the project.
Go back to Contents.
The conceptual model is includes the all the physical properties and stresses which affect the system.
- You will need to have stratigraphic data to identify boundaries between hydrologic units.
- You will need to identify the boundaries of the system. These may include hydrologic divides, impermeable beds and (some) discharge areas.
- You will need to know the hydrologic properties of each hydrologic unit. These include
- hydraulic conductivity or transmissivity
- specific storage or storage coefficient of confined units (in transient models)
- specific yield of unconfined units (in transient models)
- You will need to know the source of water for your model and its amount. This may be leakage from another aquifer or aquifer recharge.
- You need to identify all the locations and processes whereby water can leave the system. This includes pumpage from wells, discharge to streams or drainage ditches, evapotranspiraton and perhaps other processes.
- You need to make estimates of the rates of recharge or discharge by all processes.
- In transient simulations, you will need to know the initial conditions. Often the
best course is to start the simulation at a time when the potentiometric surface is close to either a static or steady state condition.
All of these will be entered into a numerical model either directly or indirectly.
Go back to the top of this page.
In MODFLOW, you must use grid lines to break up the study area into a series of blocks. When viewed from above, the blocks are rectangular but the boundaries between layers may follow geologic boundaries. The main packages used to define the grid are the Basic Package and one of the three possible Block-Centered Flow Packages. The grid cells need not all be the same size. However, no row or column should be more than 50% larger than either of it's neighbors. Any difference in size contributes to model error and if the discrepancy in the sizes of adjacent rows or columns is too large, those errors can become significant (Anderson and Woessner, 1992). MODFLOW can take a long time to run a model especially large, transient models. To minimize the time you need to wait, it's a good idea to use small grid cells only in those parts of the model where you need detailed results or where gradients are likely to be steep. You can expect to find steeper gradients near pumping wells and other discharge points. If the aquifer being modeled is anisotropic, the grid must be aligned with the directions of maximum and minimum hydraulic conductivity. Otherwise, you should align the grid in whatever way is most convenient. You can specify that some of the cells in the grid are inactive. They will be excluded from the model. The boundaries of the flow system must be approximated by a jagged edge. If recharge is not uniform, the areal extent of each zone of differing recharge should be kept close to their actual extents so that the amount of recharge is approximately correct.
Many of the cell properties are entered by using the array entry utilities - U2DREL , U1DREL, and U2DINT so you should become familiar with them especially U2DREL.
Go back to the top of this page.
You must use a consistent set of units throughout the model. If your grid scale is measured in km, and your time steps are measured in seconds, your heads must also be in km and your discharges must be in cubic km per second. Your results will be in error if you don't follow this rule.
Each cell in the model must be assigned hydrologic properties in one of the Block-Centered Flow Packages. The main source of data for these properties is well tests or other similar tests.
In some cases, you may also wish to use the Horizontal-Flow Barrier Package. It is used to simulate narrow barriers to horizontal flow such as faults and intrusions. When you have large grid cells and narrow flow barriers, you can use the Horizontal-Flow Barrier Package to reduce flow between adjacent grid cells without assigning a low transmissivity or hydraulic conductivity for either cell. The alternative would be to use much smaller grid cells but this would greatly increase the time required to run a model.
If compaction of fine grained interbeds is a significant source of water, you may wish to use the Interbed Storage Package (IBS) to determine the amount of compaction and its effects on heads.
Go back to the top of this page.
In transient models you need to choose an appropriate time step. According to Anderson and Woessner (1992), a good first estimate of the length of the time step in a homogeneous aquifer with square grid cells and one layer is
t = Sa^2/4T
where
- t=the approximate length of the longest time step that will not cause large errors.
- S=storativity (specific yield for unconfined aquifers, specific storage time aquifer thickness for confined aquifers
- a=the length of the side of one grid cell
- T=transmissivity
You should also try a shorter time step to see if it makes a difference in the computed heads or flows. Usually you want the time step to be as long as possible without adversely affecting accuracy.
Go back to the top of this page.
There are several types of boundaries possible in MODFLOW;
Go back to the top of this page.
No-flow boundaries are easy to implement. The boundaries of the grid are all no-flow boundaries. In addition you can specify cells as inactive in the Basic Package. The boundary between an active and an inactive cell is a no-flow boundary. (No flow boundaries are a special case of specified flow boundaries in which the flow is specified to be 0.)
Go back to Boundary Conditions
Go back to the top of this page.
There are three ways of implementing specified head boundaries. If you wish the head in a cell to remain constant throughout the simulation, use the IBOUND array in the Basic Package to specify it as a constant head cell. Then whatever head you specify as the starting head in the Basic Package will be retained throughout the simulation.
Sometimes, however, you may wish to specify a head for an individual cell but you want the head to change with time. You can use the Time-Variant Specified-Head Package for that. You can also use the General-Head Boundary Package, although for specified head boundaries, the General-Head Boundary Package is usually less satisfactory than the Time-Variant Specified-Head Package.
Go back to Boundary Conditions
Go back to the top of this page.
Specified flow boundaries can be implemented in several ways. For recharge across the water table, you can use the Recharge package. MODFLOW multiplies a user-specified recharge rate with units of length/time by the area of each cell (row width times column width) for which recharge is specified to determine the volumetric recharge rate. The Recharge package only supplies water to at most one cell in each vertical column. If you need a specified flow to two or more cells in one vertical column you can use the Well Package. In it, the user specifies the volumetric flow rate. If the flow rate you have is in units of length/time, you will have to figure out what the correct volumetric rate by multiplying your flow rate by the area of the side of the cell through which flow will occur.
Go back to Boundary Conditions
Go back to the top of this page.
The Drain, River, Stream, Evapotranspiration, and General-Head Boundary Packages can all be used for head dependent flow boundaries. The Drain, River, and Stream Packages are all similar in that flow can occur out of an aquifer into a surface water body. Water can also infiltrate into the aquifer in the River and Stream Packages but not in the Drain Package. The head in the River is constant whereas the head in the Stream can vary and the Stream can go completely dry. In both the River and Stream Packages, the rate of infiltration into the aquifer changes to a constant rate when the head in the aquifer drops below the bottom of the streambed or riverbed material. The main difference between the River package and the General-Head Boundary Package is that infiltration into the aquifer is unlimited in the General-Head Boundary Package.
In the Evapotranspiration Package, discharge from the aquifer occurs at some user-specified maximum rate when the water level is above some specific head and drops to 0 as the head drops to some other head (the "evapotranspiration extinction depth"). Setting the evapotranspiration extinction depth below the bottom of the cell in which evapotranspiration is occurring can cause Modflow to fail to converge to a solution.
Go back to Boundary Conditions
Go back to the top of this page.
You must select a solver to solve the finite difference equations. If you are using the wetting capability of BCF2 or BCF3, you should probably use the Preconditioned Conjugate Gradient solver or possibly the Strongly Implicit Procedure. If you are not using the wetting capability and you have a two dimensional cross-sectional (profile) model, the Slice-Successive Overrelaxation Package may have some advantages if you construct the model along a row rather than along a column. Otherwise, it's not clear which method is best.
Go back to the top of this page.
Once you have prepared the input files you need to give each file a meaningful name. Usually I use the root to identify the model run and the extension to identify the package with which the file is associated. Some people use the opposite convention. If you are using a commercial version of MODFLOW you will probably be prompted for each of the files by the name of the package to which it belongs or, in the case of output files, by their function. If you are using one of the free versions of MODFLOW, you may be prompted for some of the file names by unit number. You assign unit numbers for input files in the Basic package. The unit numbers for the output files are specified in the Output Control option and near the beginnings of most of the other packages. Keep a list available of the unit number you assigned to each input and output file. It's easiest to run the model by putting all the input files in their own directory and starting the program from there. You can have all the output files go to that directory too. The output files can be voluminous. In one recent case, the output from one of my model runs would take up approximately 4Mb. With >50 runs of similar models it was easy to fill up my hard drive. You can use data compression programs such as PKZIP or WinZip to compress the files. In some cases even the compressed file may be too large to fit on a single floppy disk. Some data compression programs such as PKZIP can split a large file among several disks. WINZIP can do this but only by acting as a shell for PKZIP.
Once you have run the model, check the water budget for each time period. Water budget errors greater than 1% on any time step are a sign of trouble. Try using a more stringent convergence criterion.
You can now compare your results with observations to see how well they compare. It's important to remember when interpreting MODFLOW output files, that the heads given for cells containing pumping wells are not the heads at those wells. The modeled heads more closely resemble the average head in the cell. Typically, the head observed at a pumping well is an extreme head; either a minimum head for pumping wells or a maximum head for recharge wells. Anderson and Woessner (1992) give the following formulas (Eq. 5.1 and 5.7) relating head in the well to modeled head in its cell for confined and unconfined aquifers.
For confined aquifers, hw = hm - (Q/2{pi}T)*ln(re/rw).
For unconfined aquifers, hw = SQRT[hm^2-(Q/{pi}K)*ln(re/rw)].
Where
- hw = head in the well,
- hm = head in the model cell,
- Q = pumping rate,
- T = aquifer transmissivity,
- re = the effective well block radius,
- rw = the well radius, and
- K = hydraulic conductivity.
For models which have a regular grid in the vicinity of the well re=0.208a where a=the length of the side of the cell. For other geometries, see Anderson and Woessner (1992).
You should keep careful records of your modeling effort as outlined by Anderson and Woessner (1992).
Go back to the top of this page.
After comparing your model to observations you will probably find some discrepancies between the model and the observations. These may be the result of observational error, inaccuracies in the conceptual model, or failure of the model to converge. The wetting capability of BCF2 and BCF3 can cause problems with convergence. Even if the model finds a solution, you should check to make sure that no cells are inactive at unreasonable times. Inactive cells act as no-flow boundaries. Such no-flow boundaries may occur at places or times where they shouldn't if cells go dry unexpectedly.
Errors in the conceptual model are probably the biggest source of error. Models are calibrated by changing input parameters such as hydraulic conductivity, transmissivity, or storativity. Errors in these values are important but there may also be errors in the locations of boundaries and incorrectly specified boundary parameters or unrecognized hydrologic processes. In attempting to calibrate a model, all aspects of the model must be considered uncertain. Even if the model matches the observations, that doesn't mean that the model is correct. It may correctly predict results for one set of input conditions but not for another. A model that works well for normal amounts of recharge for example might give incorrect results when recharge was low or high.
Finally some error can be due to observational error. The error in the measurement of the distance from ground surface or well top to the water level in the well is likely to be small. Sometimes, however, the ground elevations are not known precisely. If ground elevations are only known to within 5 m, then a 3 m error in the modeled head is insignificant. Differences in heads within a set of nested piezometers will probably be known with much less uncertainty and may make a valuable calibration target. Baseflow in streams can also be a valuable calibration target.
Go back to the top of this page.
Go back to Contents.
©1995,1996, Richard B. Winston,
rwinsto@lsuvm.sncc.lsu.edu
Dept. of Geology and Geophysics
Louisiana State University
Baton Rouge, La 70803
U.S.A.
504-388-2337
Fax 504-388-2302
http://scribe.geol.lsu.edu/rbwinston.html
ftp://aapg.geol.lsu.edu/pub/winston